Reconstruction of silicon surfaces: a stochastic optimization problem 
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Over the last two decades, scanning tunnelling microscopy (STM) has become one of the most 
important ways to investigate the structure of crystal surfaces. STM has helped achieve remarkable 
successes in surface science such as finding the atomic structure of Si(lll) and Si(OOl). For high- 
index Si surfaces the information about the local density of states obtained by scanning does not 
translate directly into knowledge about the positions of atoms at the surface. A commonly accepted 
strategy for identifying the atomic structure is to propose several possible models and analyze their 
corresponding simulated STM images for a match with the experimental ones. However, the number 
of good candidates for the lowest-energy structure is very large for high-index surfaces, and heuristic 
approaches are not likely to cover all the relevant structural models. In this article, we take the 
view that finding the atomic structure of a surface is a problem of stochastic optimization, and we 
address it as such. We design a general technique for predicting the reconstruction of silicon surfaces 
with arbitrary orientation, which is based on parallel-tempering Monte Carlo simulations combined 
with an exponential cooling. The advantages of the method are illustrated using the Si(105) surface 
as example, with two main results: (a) the correct single-step rebonded structure [e.g., Fujikawa 
et ai, Phys. Rev. Lett. 88, 176101 (2002)] is obtained even when starting from the paired-dimer 
model [Mo et al, Phys. Rev. Lett. 65, 1020 (1990)] that was assumed to be correct for many years, 
and (b) we have found several double-step reconstructions that have lower surface energies than any 
previously proposed double-step models. 
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PACS numbers: 68.35.Bs, 68.47.Fg, 68.18.Fg 



I. INTRODUCTION 

Silicon surfaces are the most intensely studied crys- 
tal surfaces since they constitute the foundation of the 
biUion-doUar semiconductor industry. Traditionally, the 
low-index surfaces such as Si(OOl) are the widely used 
substrates for electronic device fabrication. With the 
advent of nanotechnology, the stable high-index sur- 
faces of silicon have now become increasingly important 
for the fabrication of quantum devices at length scales 
where lithographic techniques are not applicable. Ow- 
ing to their grooved or faceted morphology, some high- 
index surfaces can be used as templates for the growth 
of self-assembled nanowires. Understanding the self- 
organization of adatoms on these surfaces, as well as 
their properties as substrates for thin film growth, re- 
quires atomic-level knowledge of the surface structure. 
Whether the surface unit cells are small [e.g., Si(113)] 
or large [such as Si(5 5 12)], in general the atomic- 
scale models that were first proposed were subsequently 
contestedi2iiiiSiSiL& the potential importance of stable 
Si surfaces with certain high-index orientations sparked 
many independent investigations, which led to different 
proposals in terms of surface structure. 

One of the most puzzling cases has been the (105) sur- 
face, which appears on the side-facets of the pyramidal 
quantum dots obtained in the strained layer epitaxy of 
Ge or Sii_^Ge^ {x > 0.2) on Si(OOl). Using STM imag- 
ing, Mo and coworkers proposed the first model for this 
surfaceji which was based on unrebonded monatomic 
steps separated by small (two-dimer wide) Si(001)-2xl 



terraces. Subsequently, Khor and Das Sarma reported 
another possible (105) structure with a lower density 
of dangling bonds. ^ However, the relative surface en- 
ergy of the two different reconstructions^^ was not com- 
puted, and the structure proposed in Ref. |^ had not, 
at the time, replaced the widely accepted modeli of Mo 
et al. Only very recently it has been show n^d°d^ '^^'^^ 
that the actual (105) structure is made of single-height 
rebonded steps (SR), which are strongly stabilized by 
the compressive strains present in the Ge films deposited 
on Si(001>iiiiS or Si(105)i^i^°i" Other high-index sur- 
faces such as Si(113) and Si(5 5 12) have sagas of their 
own^i^i^iii^ and only in the former case there is now 
consensus^ about the atomic structure. 

The difficulty of finding the atomic structure of a sur- 
face is not related to the resolution of the STM tech- 
niques, or to understanding of the images obtained. After 
all, it is well-known that STM gives information about 
the local density of states at the surfaces and not nec- 
essarily about atomic coordinatesi^ A common proce- 
dure for finding the reconstructions of silicon surfaces 
consists in a combination of STM imaging and electronic 
structure calculations as follows. Starting from the bulk 
truncated surface and taking cues from the experimental 
data, one proposes several atomic models for the surface 
reconstructions. The proposed models are then relaxed 
using density-functional or tight-binding methods, and 
STM images are simulated in each case. At the end of 
the relaxation, the surface energies of the structural mod- 
els are also calculated. A match with the experimental 
STM data is identified based on the relaxed lowest-energy 
structures and their simulated STM images. This proce- 
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dure has long become standard and has been used for 
many high-index orientationsi2i2iii2ii2iiLi4 As described, 
the procedure is heuristic, since one needs to rely heavily 
on physical intuition when proposing good candidates for 
the lowest energy structures. In the case of stable high- 
index Si surfaces, the number of possible good candidates 
is rather large, and may not be exhausted heuristically; 
thus, worst-case scenarios in which the most stable mod- 
els are not included in the set of "good candidates" are 
very likely. On one hand, it has been recognized^ that 
the minimization of surface energy for semiconductor sur- 
faces is not controlled solely by the reduction of the dan- 
gling bond density, but also by the amount of surface 
stress caused in the process. On the other hand, intu- 
itive reasoning can tackle (at best) the problem of lower- 
ing the number of dangling bonds, but cannot account for 
the increase in surface stress or for the possible nanoscale 
faceting of certain surfaces. For this reason, we adopt 
the view that finding the structure of high-index Si sur- 
faces is a problem of stochastic optimization, in which the 
competition between the saturation of surface bonds and 
the increase in surface stress is intrinsically considered. 

To our knowledge, a truly general and robust way of 
predicting the atomic structure of semiconductor surfaces 
-understood as finding the atomic configuration of a sur- 
face of any arbitrary crystallographic orientation without 
experimental input, has not been reported. It is not clear 
that such robust atomic-scale predictions about semicon- 
ductor surfaces can even be ventured, since theoretical 
efforts have been hampered by the lack of empirical or 
semiempirical potentials that are both fast and transfer- 
able for surface calculations. However, the long process 
which lead to the discovery of the reconstruction of the 
(105) surfaceSiSiSiiSiiiiiSii^ indicates a clear need for a 
search methodology that does not rely on human intu- 
ition. The goal of this article is to present a strategy 
for finding the lowest-energy reconstructions for an ele- 
mental crystal surface. While we hope that this strategy 
will become a useful tool for many surface scientists, the 
extent of its applicability remains to be explored. Our 
initial efforts will be focused on the surfaces of silicon 
because of their utmost fundamental and technological 
importance; nonetheless, the same strategy could be ap- 
plied for any other material surfaces provided suitable 
models for atomic interactions are available. 



II. THE MONTE CARLO METHOD 

A. General considerations 

In choosing a methodology that can help predict the 
surface reconstructions, we have taken into account the 
following considerations. First, the number of atoms 
in the simulation slab is large because it includes sev- 
eral subsurface layers in addition to the surface ones. 
Moreover, the number of local minima of the poten- 
tial energy surface is also large, as it scales roughly 




FIG. 1: Schematic computational cell: the "hot" atoms (gray) 
are allowed to move, while the bottom ones (black) are kept 
fixed at their bulk locations. Different maximum displace- 
ments As and At are allowed for the atoms that are closer to 
the surface and deeper in the bulk, respectively. 



exponentialljiii^ with the number of atoms involved in 
the reconstruction; by itself, such scaling requires the 
use of fast stochastic search methods. Secondly, the 
calculation of interatomic forces is very expensive, so 
the method should be based on Monte-Carlo algorithms 
rather than molecular dynamics. Lastly, methods that 
are based on the modification of the potential energy 
surface (PES) (such as the basin-hoping algorithm^^), al- 
though very powerful in predicting global minima, have 
been avoided as our future studies are aimed at predict- 
ing not only the correct lowest-energy reconstructions, 
but also the thermodynamics of the surface. These con- 
siderations prompted us to choose the parallel-tempering 
Monte Carlo (PTMC) algorithm22i2i for this study. Be- 
fore describing in detail the computational procedure and 
its advantages, we pause to discuss the computational cell 
and the empirical potential used. 

The simulation cell has a single-face slab geometry 
with periodic boundary conditions applied in the plane 
of the surface (denoted xy), and no periodicity in the di- 
rection (z) normal to the surface (refer to Fig. ^ . The 
"hot" atoms from the top part of the slab (corresponding 
to a thickness of 10-15 A) are allowed to move, while the 
bottom layers of atoms are kept fixed to simulate the un- 
derlying bulk crystal. Though highly unlikely during the 
finite time of the simulation performed, the evaporation 
of atoms is prevented by using a wall of infinite energy 
that is parallel to the surface and situated 10 A above it; 
an identical wall is placed at the level of the lowest fixed 
atoms to prevent the (theoretically possible) diffusion of 
the hot atoms through the bottom of the slab. The area 
of the simulation cell in the xy-plane and the number 
of atoms in the cell are kept fixed during each simula- 
tion; as we shall discuss in section IV, these assumptions 
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are not restrictive as long as we consider all the relevant 
values of the number of atoms per area. Under these 
conditions, the problem of finding the most stable recon- 
struction reduces to the global minimization of the total 
potential energy V (x) of the atoms in the simulation cell 
(here x denotes the set of atomic positions). In order to 
sort through the numerous local minima of the poten- 
tial V{x), a stochastic search is necessary. The general 
strategy of such search (as illustrated, for example, by 
the simulated annealing technique^SiS^) is to sample the 
canonical Boltzmann distribution exp [— T^(x) / {ksT)] for 
decreasing values of the temperature T and look for the 
low-energy configurations that are generated. 

In terms of atomic interactions, we are constrained to 
use empirical potentials because the highly accurate ab- 
initio or tight-binding methods are prohibitive. Since 
this work is aimed at finding the lowest energy recon- 
structions for arbitrary surfaces, the choice of the empir- 
ical potential is crucial, as different interaction models 
can give different energetic ordering of the possible re- 
constructions. Furthermore, the true structure of the 
surface may not even be a local minimum of the poten- 
tial chosen to describe the interactions: it is the case, 
for example, of the adatom-interstitial reconstructions^ 
of Si(113), which are not local minima of the Stillinger- 
Weber potential.^'' The work of Nurminen et alm^ in- 
dicates that the most popular empirical potentials for 
siliconi2iiS& are not suitable for finite-temperature simu- 
lations of surfaces. After thorough numerical experimen- 
tation with several empirical potentials, we have chosen 
to use the highly optimized empirical potential (HOEP) 
recently developed by Lenosky et al?'^ HOEP is fitted to 
a database of ab-initio calculations that includes struc- 
tural and energetic information about small Si clusters, 
which leads to a superior transferability to the different 
bonding environments present at the surfaced 



B. Advantages of the parallel tempering algorithm 
as a global optimizer 

The parallel tempering Monte Carlo method (also 
known as the replica-exchange Monte-Carlo method) 
consists in running parallel canonical simulations of many 
statistically independent replicas of the system, each at 
a different temperature Ti < T2 < . . . < Tjy. The set 
of N temperatures {Ti, i = 1,2,...N} is called a tem- 
perature schedule, or schedule for short. The probability 
distributions of the individual replicas are sampled with 
the Metropolis algorithm,^* although any other ergodic 
strategy can be utilized. The key feature of the paral- 
lel tempering method is that swaps between replicas of 
neighboring temperatures Ti and Tj (j = i ± 1) are pro- 
posed and allowed with the conditional probabilitjiSSiSi 
given by 



where V{xi) represents the energy of the replica i and ks 
is the Boltzmann constant. The conditional probability 
^ ensures that the detailed balance condition is satisfied 
and that the equilibrium distributions are the Boltzmann 
ones for each temperature. 

In the standard Metropolis samplingSS of Boltzmann 
distributions, the probability that the Monte Carlo 
walker escapes from a given local minimum decreases 
exponentially as the temperature is lowered. In turn, 
the average number of Monte Carlo steps needed for the 
walker to escape from the trapping local minimum in- 
creases exponentially with the decrease of the tempera- 
ture, a scaling that makes the search for a global mini- 
mum inefficient at low temperatures. To cope with this 
problem, the parallel tempering algorithm takes advan- 
tage of the fact that the Metropolis walkers running at 
higher temperatures have larger probabilities of jump- 
ing over energy barriers. Parallel tempering significantly 
decreases the time taken for the walker to escape from 
local minima by providing an additional mechanism for 
jumping between basins, namely the swapping of configu- 
rations between replicas running at neighboring temper- 
atures. Therefore, if a given (low-temperature) replica of 
the system is stuck in a local minimum, the configura- 
tion swaps with walkers at higher temperatures can pro- 
vide that replica with states associated with other basins 
(wells on the potential energy surface) , ultimately driving 
it into the global minimum. 

Because of this swapping mechanism, parallel temper- 
ing enjoys certain advantages (as a global optimizer) 
over the more popular simulated annealing algorithm 
(SA)i22iS In order for SA to be convergent (i.e. to reach 
the global optimum as the temperature is lowered) the 
cooling schedule must be of the form22i2£ 

^^^w^Tm' (2) 

log(i + to) 

where Tq and io are sufficiently large constants. Such a 
logarithmic schedule is too slow for practical applications, 
and faster schedules are routinely utilized. Common SA 
cooling schedules, such as the geometric or the linear 
onesr^ make SA non-convergent: the Monte Carlo walker 
has a non-zero probability of getting trapped into minima 
other than the global one. 

The cooling schedule implied by Eq. |2l is, of course, 
asymptotically valid in the limit of low temperatures. In 
the same limit, the PT algorithm allows for a geometric 
temperature schedule^ii^S When the temperature drops 
to zero, the system is well approximated by a multidimen- 
sional harmonic oscillator and the acceptance probability 
for swaps attempted between two replicas with temper- 
atures T < T' is given by the incomplete beta function 
law2^ 

Ac{T,T') ~ -^rrnr-rn^ / e-^/^-^il-ef/^-^de , 



/T,-l/T.)[V(x,)-y(x,)]/fc£ 



(1) 



Bid/2, d/2) Jo 

(3) 

where d denotes the number of degrees of freedom of the 
system, B is the Euler beta function, and R = T' /T. 
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FIG. 2: Heat capacity of a Si(105) slab plotted as a function of 
temperature. The peak is located at 1550K; in order to avoid 
recalculation of the heat capacity for systems with different 
numbers of atoms and surface orientations, we set Tmax ~ 
1600K as the upper limit of the temperatures range used in 
the PTMC simulations. 

Since it depends only on the temperature ratio R, the 
acceptance probability (PJ has the same value for any 
arbitrary replica running at a temperature Ti, provided 
that its neighboring upper temperature Ti+i is given by 
Ti+i = RTi. The value of R is determined such that the 
acceptance probability given by Eq. © attains a pre- 
scribed value p, usually chosen greater that 0.5. Thus, 
the (optimal) schedule that ensures a constant probabil- 
ity p for swaps between neighboring temperatures is a 
geometric progression: 

T, = R'-'Trmn, l<i<N, (4) 

where Tmin = Ti is the minimum temperature of the 
schedule. Though more research is required in order to 
better quantify the relative efficiency of the two differ- 
ent algorithms SA and PT, it is apparent from Eqs. Q 
and Q that the parallel tempering algorithm is a global 
optimizer superior to SA because it allows for a faster 
cooling schedule. Direct numerical comparisons of the 
two methods have confirmed that parallel tempering is 
the superior optimization technique^ The ideas of par- 
allel tempering and simulated annealing are not mutually 
exclusive, and in fact they can be used together to de- 
sign even more efficient stochastic optimizers. As shown 
below, such a strategy that combines parallel tempering 
and simulated annealing is employed for the present sim- 
ulations. 



C. Description of the algorithm 

The typical Monte Carlo simulation done in this work 
consists of two main parts that are equal in terms of com- 
putational effort. In the first stage of the computation. 



we perform a parallel tempering run for a range of tem- 
peratures [T^in, Tmax]- The configurations of minimum 
energy are retained for each replica, and used as start- 
ing configurations for the second part of the simulation, 
in which each replica is cooled down exponentially until 
the largest temperature drops below a prescribed value. 
As a key feature of the procedure, the parallel temper- 
ing swaps are not turned off during the cooling stage. 
Thus, we are using a combination of parallel tempering 
and simulated annealing, rather than a simple cooling. 
Finally, the annealed replicas are relaxed to the nearest 
minima using a conjugate-gradient algorithm. We now 
describe in detail the stochastic minimization procedure. 
We shall focus, in turn, on discussing the Monte Carlo 
moves, the choice of the temperature range \Tmim Tmax], 
and the total number of replicas iV. 

The moves of the hot atoms consist in small random 
displacements with the x, y, z components given by 

A{2Ua - 1) 

where Ua {a — x,y,z) are independent random 
variables^ uniformly distributed in the interval [0,1], 
and A is the maximum absolute value of the displace- 
ment. We update the positions of the individual hot 
atoms one at a time in a cyclic fashion. Each attempted 
move is accepted or rejected according to the Metropolis 
logic. A complete cycle consisting in attempted moves 
for all hot particles is called a pass (or sweep) and con- 
stitutes the basic computational unit in this work. We 
have computed distinct acceptance probabilities for the 
hot atoms that are closer to the surface (situated within 
a distance of 5 A below the surface) and for the deeper 
atoms, the movements of which are essentially small os- 
cillations around the equilibrium bulk positions. Conse- 
quently, as shown in Fig. ^ we have employed two dif- 
ferent maximal displacements, for the surface atoms, 
and Af, for the bulk-like atoms lying in the deeper sub- 
surface layers. The displacements As and Af, are tuned 
in the equilibration phase of the simulation in such a way 
that the Monte Carlo moves are accepted with a rate of 
40% to 60%. This tuning of the maximal displacements 
has been performed automatically by dividing the equili- 
bration phase into several blocks, computing acceptance 
probabilities for each block, and increasing or decreasing 
the size of the displacements A^^f, until the acceptance 
probabilities reached values between 40% and 60%. The 
automatization is necessary because the optimal displace- 
ments computed for replicas running at different temper- 
atures have different values. The maximal displacement 
As for the surface atoms is found to be larger than the 
maximal displacement for the bulk-like atoms. Though 
expected in view of the larger mobility of the surface 
atoms, the difference between A^ and Af, is not substan- 
tial and the reader may safely employ a single maximal 
displacement for all hot atoms at a given temperature. 

Parallel tempering configuration swaps are attempted 
between replicas running at neighboring temperatures at 
every 10 passes in an alternating manner, first with the 
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closest lower temperature then with the closest higher 
temperature. Exception make the two replicas that run 
at end temperatures Ti = Tmin and Tn = Tmax, which 
are involved in swaps every 20 passes. The range of 
temperatures [Tmin,Tmax] and the temperature sched- 
ule Ti < T2 < • • • < T/v have been chosen as described 
below. 

The maximum temperature Tmax must be high enough 
to ensure that the corresponding random walker has good 
probability of escaping from various local minima. How- 
ever, as the temperature is raised, increasingly more ther- 
modynamic weight is placed on local minima that have 
high energies compared to the global minimum. Stillinger 
and Weberi?iiS have argued that the number of local 
minima increases exponentially with the dimensionality 
of the system. As such, the probability that the walker 
visits the basin of the global minimum significantly de- 
creases with the increase of temperature. A very strong 
decrease occurs at the melting point, beyond which most 
of the configurations visited are associated with the liq- 
uid phase. The basins of these configurations are un- 
likely to contain the global minimum or, in fact, any of 
the low-energy local minima associated with meaningful 
surface reconstructions. Therefore, the high-temperature 
end must be set equal to the melting temperature. 

The melting temperature of the surface slab can be de- 
termined from a separate parallel tempering simulation 
by identifying the peak of the heat capacity plotted as 
a function of temperature. As Fig. El shows, the melting 
temperature of a Si(105) sample slab with 70 hot atoms 
is about 1550 K. Rather than determining a melting 
temperature for each individual system studied, we have 
employed a fixed value of Tmax = 1600 K. The melt- 
ing temperature of the slab determined here (Fig. |2l is 
different from the value of 1250K reported for the bulk 
crystali^ the discrepancy is due to surface effects, finite- 
size effects, as well as to the fact that the hot atoms are 
always in contact with the rigid atoms from the bottom 
of the slab. Though we use Tmax = 1600K for all sim- 
ulations, we note that differences of 100K-200K in the 
melting temperature of the slab do not significantly af- 
fect the quality of the Monte Carlo sampling. For most 
surfaces and system sizes of practical importance, the 
value of 1600 K is in fact un upper bound for the melting 
temperature; this may sometimes cause the one or two 
walkers that run at the highest temperatures to be un- 
coupled from the rest of the simulation, since they might 
sample amorphous or liquid states. However, this loss in 
computational resources is very small compared to the 
additional effort that would be required by a separate 
determination of the heat capacity for each surface slab 
used. 

In theory, the lowest temperature Tmin should be set 
so low that the walker associated with this temperature is 
virtually localized in the basin associated with the global 
minima. Nevertheless, obstacles concerning the efficient 
use of computational resources prevent us from doing so. 
Numerical experimentation has shown that a tempera- 
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FIG. 3: Exponential cooling of the A'^ = 32 Monte Carlo walk- 
ers (replicas of the surface slab) used in the simulation. For 
clarity, only eight walkers are shown (every fourth walker). 
The cooling is performed in 18 steps: at each step the tem- 
perature is modified by the same factor a = 0.85 for all walk- 
ers, Eq. ©. For every cooling step k, we have a different 
parallel tempering schedule where each replica is coupled to 
the walkers running at neighboring temperatures via config- 
uration swaps [Eq. I^J with R = 4^'''^^]. This coupling is 
symbolized by the double-arrow lines in the inset. 



ture of Tmin = 400 K is low enough that only local min- 
ima associated with realistic surface reconstructions are 
frequently visited. A further selection among these local 
minima is performed in the second part of the Monte 
Carlo simulation, when all temperatures of the initial 
schedule {Ti, i — 1,2, ...N} are gradually lowered to val- 
ues below 100 K; as it turns out, this combination of 
parallel tempering and simulated annealing makes opti- 
mal use of computational resources. Below the melting 
point the heat capacity of the surface slab is almost con- 
stant and well approximated by the capacity of a multidi- 
mensional harmonic oscillator (refer to Fig. [3). In these 
conditions, the acceptance probability for swaps between 
neighboring temperatures T and T' is given by Eq. (O 
(see also Ref. 132) . It follows that the optimal tempera- 
ture schedule on the interval [Tmin , Tmax] is the geometric 
progression Q), where 



R — {Tmax /Tmin) 



l/[N{d,p)-l] 



We have written N = N{d,p) to denote the smallest 
number of replicas that guarantees a swap acceptance 
probability of at least p for a system with d degrees 
of freedom. Since the best way to run PTMC calcula- 
tions is to use one processor for each replica of the sys- 
tem, the feasibility of our simulations hinges on values of 
N{d,p) that translate directly into available processors. 
The number of walkers N{d,p) can be estimated^ by 



N{d,p) 



4crf"^(l-p) 



(5) 
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where [x] denotes the largest integer smaller than x, and 
erf~^ is the inverse error function. Based on Eq. we 
have used = 32 walkers for all simulations, which en- 
sures a swap acceptance ratio greater than p — 0.5 for 
any system with less than 300 hot atoms, d < 900. The 
first part of all Monte Carlo simulations performed in the 
present article consists of a number of 36 x 10^ passes for 
each replica, preceded by 9 x 10'* passes allowed for the 
equilibration phase. When we retained the configura- 
tions of minimum energy, the equilibration passes have 
been discarded so that any memory of the starting con- 
figuration is lost. 

We now describe the second part of the Monte Carlo 
simulation, which consists of a combination of simu- 
lated annealing and parallel tempering. At the k-th 
cooling step, each temperature from the initial temper- 
ature schedule {Ti,i — 1,2,..N} is decreased by a fac- 
tor which is independent of the index i of the replica, 
Tj^*^-* = akT'^^~^\ Because the parallel tempering swaps 
are not turned off, we require that at any cooling step k 
all N temperatures must be modified by the same fac- 
tor Uk in order to preserve the original swap acceptance 
probabilities. The specific way in which depends on 
the cooling step index k is determined by the kind of an- 
nealing being sought. In this work we have used a cooling 
schedule of the form 

where T, = T^^^ and a is determined such that the 
temperature intervals [T^''~'^\t^^~^'^] and [rf ^\ T^''^] 
spanned by the parallel tempering schedules before and 
after the k-th cooling step overlap by 80%. This yields a 
value for a given by (0.2T„,„ + O.STmax) /Tmax = 0.85. 
We have also tested values of a larger than 0.85, and did 
not find any significant improvements in the quality of 
the sampling. 

The reader may argue that the use of an exponential 
annealing [Eq. ©] is not the best option for attaining 
the global energy minimum of the system. Apart from 
the theoretical considerations discussed in the preceding 
subsection that only a logarithmic cooling schedule would 
ensure convergence to the ground state^SiS&it is known 
that the best annealing schedules for a given computa- 
tional effort oftentimes involve several cooling-heating cy- 
cles. We emphasize that in the present simulations, the 
most difficult part of the sampling is taken care of by the 
initial PTMC run. In addition, since the configuration 
swaps are not turned off during cooling (refer to Fig.|2J), 
the Monte-Carlo walkers are subjected to cooling-heating 
cycles through the parallel tempering algorithm. 

The purpose of the annealing (second part of the simu- 
lation) is to cool down the best configurations determined 
by the initial parallel tempering in a way that is more ro- 
bust than the mere relaxation into the nearest local min- 
imum. If the initial PTMC run is responsible for placing 
the system in the correct funnels (groups of local minima 
separated by very large energy barriers), the annealing 



part of the simulation takes care of jumps between lo- 
cal minima separated by small barriers within a certain 
funnel. For this reason, the annealing is started from the 
configurations of minimum energy determined during the 
first part. The cooling is stopped when the largest tem- 
perature in the parallel tempering schedule drops below 
lOOK. This criterion yields a total of 18 cooling steps, 
with 2 X 10'' MC passes per replica performed at every 
such step. 

Each cooling step is preceded by 5 x 10'^ equilibration 
passes, which are also used for the calculation of new 
maximal displacements and A{,, as these displace- 
ments depend on temperature and must be recomputed. 
In fact, each cooling step is a small-scale version of the 
first part of the simulation. The only difference is that 
the cooling steps are not started from the configurations 
of minimum energy determined at the preceding cool- 
ing steps. (Otherwise, because the number of passes for 
a given step is quite small, the walkers might not have 
time to escape from some spurious local minima and we 
would end up restarting them over again from the respec- 
tive minima.) 

The third and final part of the minimization procedure 
is a conjugate-gradient optimization of the last configu- 
rations attained by each replica. The relaxation is nec- 
essary because we aim to classify the reconstructions in 
a way that does not depend on temperature, so we com- 
pute the surface energy at zero Kelvin for the relaxed 
slabs i, j = l,2,...7V. The surface energy 7 is defined as 
the excess energy (with respect to the ideal bulk config- 
uration) introduced by the presence of the surface: 



where Em is the potential energy of the atoms that 
are allowed to move, Cf, = 4.6124eV is the bulk cohesion 
energy given by HOEP, and A is the surface area of the 
slab. 



III. RESULTS FOR THE SI(105) SURFACE 

We have tested the method for a variety of surface 
orientations, such as (113), (105) and (5 5 12). In this 
section we are presenting results for Si(105), a choice that 
was determined by the ubiquity of the (105) orientation 
on the side facets of the pyramidal quantum dots ob- 
tained in the heteroepitaxial deposition of Ge and Si-Ge 
alloys on Si(OOl). Recent experimental and theoretical 
work on the atomic structure of (105) surfaces^iifliiLiSii^ 
provides a strong testing ground for the current investiga- 
tions. In order to assess the versatility of the method and 
to provide a direct comparison with a previous heuristic 
studjii^ of the (105) reconstructions, we start our PTMC 
simulations from each of the structures found in Ref . [3- 
To establish the nomenclature for the discussion to fol- 
low, we recall that the structures were labelled by SU, 
SR, DU, DUl, DR, DRl, and DR2, where the first letter 
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FIG. 4: Total energies of the 32 replicas of the Si(105) computational slab at the end of the cooling sequence (circles), and after 
the subsequent conjugate-gradient relaxation (triangles). The PTMC procedure has been started with all the replicas in the 
same configuration taken from the set reported in Ref. [ij: SU(a), DU(b), DR(c), DR2(d). The lowest-energy configurations 
depend on the total number of atoms n, which is indicated in each panel. Six new double-step structural models are found, 
denoted by DT, DXl, DX2, DR2a, DR2/9 and DR27, with surface energies smaller than those of the corresponding starting 
structures 




SU DT SR 



FIG. 5; Si(105) reconstructions obtained when starting from the SU model: SU, DT and SR. The DT structure is a novel 
double-stepped structure retrieved by replicas running at intermediate temperatures (see also Fig. BI^))- The single-step 
rebonded structure2*Si— (SR) is the global optimum. The rectangle represents the surface unit cell, which is the same as 
the as the periodic supercell used in the simulations. Atoms are rainbow-colored according to their coordinate along the [105] 
direction, with the red atoms being at the highest positions. 
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denotes the height of the steps (single S, or double D), the 
second letter indicates whether the step is rebonded (R) 
or unrebonded (U) and the digit distinguishes between 
different structures that have the same broad topological 
features.^"' These reconstructions have different numbers 
of atoms and different linear dimensions of the periodic 
cell. The dimensions of the cell are chosen 2a x 
(a = 5.431A is the bulk lattice constant of Si) for all the 
models considered except DR2, whose topology requires 
a periodic cell of 2a x 2a\/6.5. The thickness of the slab 
corresponds to two unit cells in the z direction, with a 
maximum of 208 atoms, of which only about half are 
allowed to move. 

The results of the PTMC simulations for the Si(105) 
surface are plotted in Fig. 2) which shows the total energy 
for each of the = 32 replicas at the end of the cool- 
ing procedure (circles) and after the conjugate-gradient 
relaxation (triangles). Figs. ^Ja), (b), (c) and (d) show 
the total energies of the reconstructions obtained starting 
from the SU, DU, DR and DR2 models, respectively. In 
each case, we have obtained at least two structures with 
lower surface energies than the starting configurations, 
which we discuss in turn. 

Fig. ^a) shows that the (starting) SU structure'*' is 
found only by the two replicas running at the highest 
temperatures, while colder walkers find a novel double- 
stepped structure, termed here "transitional" (DT). At 
even lower temperatures, the double-steps of the DT re- 
construction unbunch into single-height rebonded (SR) 
steps; the three different configurations that correspond 
to the energies plotted in Fig. ^a) are shown in Fig. |S1 
Therefore, the correct SR structure^iiSiiiiiSii^ is retrieved 
even when starting from the topologically different SU 
model. The usefulness of this PTMC simulation becomes 
apparent if we recall that the SU structure was widely be- 
lieved to be correct for more than a decade after its publi- 
cation. As we shall see, the ground state obtained in our 
stochastic search is independent of the initial configura- 
tion. The only condition for finding the reconstruction 
with the lowest surface energy is to prescribe the correct 
number of atoms and the correct dimensions for the sim- 
ulation slab. We will address these practical aspects in 
the next section; for now, we continue to describe the 
results obtained for different numbers of atoms in the 
computational slab. 

The simulation that starts from the DU model finds 
two distinct rebonded structures, denoted by DXl and 
DX2 in Fig.^b). Both these structures are characterized 
by the presence of single dimers at the location of steps 
(see Fig. IS)), which reduces the number of dangling bonds 
per unit area from 6db/a'^y/K5 (starting structure DU) to 
5db/a^^/6.5. The DXl reconstruction is the most stable, 
and it is obtained in all but three replicas of the system. 

Although it has a small density of dangling bonds, the 
DR structure has large surface energy due to the V2 x 1 
terrace reconstructionii^ Since the density of dangling 
bonds is the lowest possible (4c?&/a^\/6.5), the minimiza- 
tion of surface energy is dictated by the reduction of sur- 



face stress. Unlike the case of SU and DU structures 
(described above), not a single replica have retained the 
starting model DR. Instead, the DT and SR structures 
are retrieved (refer to Fig. Efc)). When starting from 
the DR2 structure we obtain at least three low energy 
structures denoted by DR2a, DR2/3 and DR27 ( Fig . 171) 
which have not been previously proposed in Refs. Il3.l36l 
or elsewhere. Owing to a larger area of the slab, por- 
tions of the newly reconstructed unit cells have patches 
that resemble the models obtained in prior simulations. 
In particular, the atomic scale features of the steps on 
DR2a are very similar to those of the SR structure, a 
similarity that reflects in the very small relative surface 
energy of the two models (w 1.6meV/A^ ). 

We note that the simulations described have a total 
number of atoms that is between n = 202 and n — 206 
(Fig. EI^) and (c)) per 2a'^^/K5 area. To cover all the 
possibilities for intermediate numbers of atoms, we also 
perform a simulation with n = 205; this value of n does 
not correspond to any of the models reported in Ref. 0, 
and the parallel tempering run is started from a bulk- 
truncated configuration. In this case two new structures 
are found; these structures are named DYl and DY2 
and shown in Fig. IHl [The letters X and Y appearing 
in DXl, DX2, DYl, DY2 (all denoting double-stepped 
rebonded structures. Fig. (H)) do not stand for particular 
words, they are simply intended to unambiguously la- 
bel the structures in a way that does not complicate the 
notation.] While for the DYl model the rebonding is re- 
alized via bridging bonds ^i^i in the case of DY2 we find 
unexpected topological features such as fully saturated 
surface atoms and over-coordinated bulk atoms. Even 
though these structural units (seen in the DY2 panel of 
Fig. reduce the number of dangling bonds, they also 
create high atomic-level stresses which make the DY2 re- 
construction relatively unfavorable. 

We have also performed PTMC simulations with SR, 
DRl and DUl as initial configurations, but have not ob- 
tained any other novel reconstructions. We found that 
SR and DRl are the global energy minima correspond- 
ing to 206 atoms and 203 atoms, respectively. The DUl 
structures^ (202 atoms) has lead to the same reconstruc- 
tions as the SU model (206 atoms). This result indicates 
a periodic behavior of the surface energy as a function of 
the total number of atoms, which will be discussed next. 



IV. DISCUSSION 

To further test that the lowest energy states for given 
number of atoms are independent of the initial config- 
urations, we have repeated all the calculations using 
bulk-truncated surface slabs (Fig. ^ instead of recon- 
structed ones. We have varied the number of atoms n 
in the simulation cell between 196 and 208, where the 
latter corresponds to four bulk unit cells of dimensions 
a\/K5 X a X a-\/6.5 stacked two by two in the [010] and 
[105] directions. For the cases with n < 208, we have 



FIG. 6: Double-step reconstructions of Si(105) with periodic cells (rectangles shown) of dimensions 2a x a\/6.5. The color 
scheme is the same as in Fig. El Except for DRl, all other structures are new. 



started the PTMC simulations from structmes obtained 
by taking out a prescribed number atoms from random 
surface sites, and have found the same ground state ir- 
respective of the locations of the removed atoms. For 
values of n equal to 202, 203, 204 and 206, the ground 
states (global minima) are also the same as the ones ob- 
tained from the reconstructed models DR, DRl, DU, and 
SU, respectively. Furthermore, we have tested that even 
when removing arbitrary subsurface atoms the simulation 
retrieves the same ground states without increasing the 
computational effort. This finding speaks for the quality 
of the Monte Carlo sampling and gives confidence in the 
predictive capabilities of the method described in section 
IL The lowest surface energies obtained at the end of the 
numerical procedure are shown in Fig. |^ as a function 
of the number of atoms in the simulation cell. As illus- 
trated in Fig. 121 the simulation finds the same ground 
states at periodic intervals of An = 4. At first sight, this 
is somewhat surprising given that the number of under- 
coordinated surface atoms in a bulk-truncated cell of di- 
mensions 2a X a\/6.5 is twelve (refer to Fig. (SJ. The re- 
duced periodicity of the surface energy with the number 
of atoms in the supercell is due to the underlying crystal 
structure, which lowers the number of symmetry-distinct 
global minima to only four. Thus, we have considered 
all possibilities in terms of numbers of atoms in a simu- 
lation slab of area 2a^\/6.5. The surface energies of the 
optimal reconstructions for relevant values of n, as well 
as those of some higher-energy structures, are collected 
in TableQ] As shown in the table, the global minimum of 
the surface energy of Si(105) is obtained for the single- 



height rebonded structure SR. While this finding is in 
agreement with recent reports (SiiSiiLiSii^ it is the result 
of an exhaustive search rather than a comparison between 
twc&iSiiiiH or morei^ heuristically proposed structures. 

From Table U we also note that the SR and the 
DR2a structures have surface energies that are within 
1.6meV/A^ from one another. This gap in the surface en- 
ergy of the two models (SR and DR2) is smaller than the 
expected accuracy of relative surface energies determined 
by an empirical potential. Therefore, it is very likely that 
these two reconstructions can both be present on the 
same surface under laboratory conditions. As recently 
pointed out,^'^ the coexistence of several configurations 
with different topological features but similar surface en- 
ergies gives rise to the atomically rough and disordered 
aspec t'^^1'^^ of the Si(105) surface. The surface energies 
computed using HOEP for various rebonded structures 
(Table are close to the values obtained previouslyi^ 
at the tight-binding leveli^ For the unrebonded struc- 
tures (SU and DU), the differences between the HOEP 
values and the tight-binding ones are larger: this dis- 
crepancy is caused by the inability of the HOEP interac- 
tion model to capture the tilting of the surface dimers, 
which is an important mechanism for the relaxation of 
these unrebonded configurations. Despite this shortcom- 
ing, we have found that the HOEP potential is accurate 
enough to predict the correct bonding topology of the 
global minimum reconstructions for a variety of surface 
orientations. If a comparison with experimental STM im- 
ages is desired, further geometry optimizations are nec- 





FIG. 7: Double-step reconstructions of Si(105) with periodic 
cells (rectangles shown) of 2a x 2a\/6.5. Although the starting 
structure [the DR2 model-- shown in (d)] has a reasonably low 
dangling bond density (5d6/a^\/6.5), the Monte Carlo sim- 
ulation has retrieved three more reconstructions, all having 
smaller surface energies (refer to Table H}. These novel struc- 
tures [shown in figs, (a)-(c)] are labelled by DR2q, DR2/3, 
DR27. The atoms are rainbow-colored as indicated in Fig. |3 

essary at the level of electronic structure methods: these 
calculations would have to consider different tiltings of 
the surface bonds, and in each case the simulated im- 
age is to be compared with the experimental one. Thus, 
even for surfaces where dimer tilting is important, the 
Monte Carlo simulation based on the HOEP interaction 
modelSi can still serve as a very efficient tool to find good 
candidates for the lowest energy structures. 

Two practical issues have to be addressed when us- 



FIG. 8: Atomic structure of the bulk truncated Si(105) sur- 
face, viewed from the side (a) and from the top (b). The 
rectangle of dimensions 2a x a\/KE marks the periodic cell 
used in most of the simulations, and contains two unit cells of 
the bulk-truncated surface. For clarity, only a single subsur- 
face (001) layer is shown. In this picture (unlike in Figs.|S]|S| 
and[7|l atoms are colored according the their number of dan- 
gling bonds (db) before reconstruction: red= 2db, green= Idb 
and blue= Odb. 

ing PTMC simulations for surface structure prediction. 
First consideration is related to the size of the compu- 
tational cell. If a periodic surface pattern exists, the 
lengths and directions of the surface unit vectors can be 
determined accurately through experimental means (e.g., 
STM). In those cases, the periodic lengths of the simula- 
tion slab should simply be chosen the same as the ones 
found in experiments. On the other hand, when the sur- 
face does not have two-dimensional periodicity (as it is 
the case of unstrained Si(105) surfac o^^i'^^ ), or when ex- 
perimental data is not available, one should systemati- 
cally test computational cells with periodic vectors that 
are low-integer multiples of the unit vectors of the bulk 
truncated surface; the latter unit vectors can be easily 
computed from the knowledge of crystal structure and 
surface orientation. Secondly, the number of atoms in 
the simulation cell is not a priori known, and there is no 
simple criterion to find the set of numbers that yield the 
lowest surface energy for a slab with arbitrary orienta- 
tion. Adapting the algorithm presented in section II for 
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ber of atoms in the computational cell is not new, as it 
also appears, for example, in classic algorithms for pre- 
dicting the bulk crystal structurei^ As shown above for 
the case of Si(105), a successful way to deal with this 
problem is to simply repeat the simulation for systems 
with consecutive numbers of atoms, and look for a peri- 
odic behavior of the surface energy of coldest replicas as a 
function of the number of particles in the computational 
cell. Note that if the thickness of the slab is sufficiently 
large, such periodicity of the lowest surface energy with 
respect to the number of atoms in the supercell is guar- 
anteed to exist: in the worst case, the periodicity will 
appear when an entire atomic layer has been removed 
from the simulation cell. 



FIG. 9: Surface energy of the global minimum structure plot- 
ted versus the total number of atoms n in the simulation 
slab. Even though there are twelve under-coordinated atoms 
in each bulk-truncated periodic cell (refer to Fig.|HJ, the values 
of the surface energy repeat at intervals of An = 4. The un- 
derlying bulk structure reduces the number of distinct global 
minima to four. 



n Structure Bond counting HOEP TB 
(db/a'^VKE) (meV/A^) (meV/A^) 



206 


SR 


4 


82.20 


82.78 




DT 


4 


85.12 






SU 


6 


88.35 


83.54 


205 


DYl 


5 


86.73 






DY2 


4.5 


88.61 




204 


DXl 


5 


84.90 






DX2 


5 


86.04 






DU 


6 


90.18 


84.84 


203 


DRl 


5 


86.52 


85.22 


2x203 


DR2a 


4.5 


83.77 




2x203 


DR2/3 


4.5 


84.64 




2x203 


DR27 


4.5 


86.15 




2x203 


DR2 


5 


86.34 


83.48 



TABLE I: Surface energies of different Si(105) reconstruc- 
tions, calculated using the HOEP interatomic potentiali— 
The structures are grouped according to the number of atoms 
n in the simulation cell. Atomic configurations of selected re- 
constructions are shown in Figs.|S]|n]andEI The third column 
shows the number of dangling bonds (db) per unit area, ex- 
pressed in units of a^V6.5. The last column indicates the 
tight-bindin^^(TB) values reported in Ref. ll3L 



a grand-canonical ensemble is somewhat cumbersome, as 
one would have to consider efficiently the combination 
of two different types of Monte Carlo moves: the small 
random displacements of the atoms (continuous) and the 
discrete processes of adding or removing atoms from the 
simulation slab. The problem of finding the correct num- 



V. CONCLUDING REMARKS 

In conclusion, we have developed and tested a stochas- 
tic method for predicting the atomic configuration of sil- 
icon surfaces. If suitable empirical models for atomic 
interactions are available, this method can be straight- 
forwardly applied for the determination of the structure 
of any crystallographic surface of any other material. 
Using the example of Si(105), we have shown that the 
PTMC search is superior to heuristic approaches because 
it ensures that the topology corresponding to the low- 
est surface energy is considered in the set of good pos- 
sible structural models. We have performed an exhaus- 
tive search for different numbers of atoms in the simula- 
tion cell and have found that the global minimum of the 
(105) surface energy is the single-height rcbonded model 
SR, in agreement with recent studieSf?'^"'^^'^^-i2, The ex- 
periments of Zhao et alm^ indicated that double-stepped 
structures are present on the unstrained Si(105) surface: 
our simulations indeed have found double-stepped mod- 
els with surface energies that are close to the surface 
energy of the optimal SR reconstruction. In addition, 
these double-stepped models (termed DR2q;, DR2/3, and 
DR27) are energetically more favorable than the double- 
stepped structures proposed in Refs. andlT^ 

We would like to comment on the key role played by 
the empirical potential used in the present simulations. 
A highly transferable interatomic potential is required for 
a satisfactory energetic ordering of different reconstruc- 
tions. While we would not expect any empirical potential 
to accurately reproduce the relative surface energies of all 
the reconstructions found, we can at least expect that the 
chosen potential correctly predicts the bonding topology 
for well-known surface reconstructions. In this respect, 
the HOEP modelSi proved superior to the most widely 
used interatomic potentialsi2^i2& Given this comparison, 
the results presented here would represent a validation 
of the workSi towards more transferable potentials for 
silicon. We also hope that these results would stimulate 
further developments of interatomic potentials for other 
semiconductors 

With the exception of Si(105), Si(113)3 and (likely) 
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Si(114)ji^ the atomic structure of other stable high-index 
sihcon surfaces has not been elucidated, although a sub- 
stantial body of STM images has accumulated to date4ii 
A similar situation exists for Ge surfaces as welliii The 
methodology presented in this article can be used (ei- 
ther directly or in combination with the STM images^) 
to determine the configuration of all high-index Si sur- 
faces, as long as the HOEP potential remains satisfac- 
tory for all orientations to be investigated. Furthermore, 
with certain modifications related to the implementa- 
tion of empirical potentials for systems with two atomic 
species, the PTMC method could help bring important 
advances in terms of finding the thermodynamically sta- 
ble intermixing composition of various nanostructures 
obtained by heteroepitaxial deposition of thin films on 
silicon substrates. Though such studies have recently 
been reported, only the intermixing at a given atomic 
bonding topology has been investigated. The interplay 



between reconstruction and intermixing is another chal- 
lenging and important problem that can be tackled via 
PTMC simulations. Lastly, the method presented in this 
article may also be used for studying the decomposition 
of unstable orientations into nanofacets, as well as for 
predicting the thermodynamics of surfaces in the pres- 
ence of adsorbates or applied strain. 
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